///@file SymTensor2.cpp
///@brief Implementation of the class representing a symetric second order tensor
///@author Arnaud Duval
///@version 0.0
///@date 2010/03/08


#include "SymTensor2.h"

///@brief Default constructor
SymTensor2::SymTensor2()
{
	dim = 0;
}

///@brief Constructor giving the dimension
///@param n dimension of the created SymTensor
SymTensor2::SymTensor2(int n)
{
	dim = n;
	a = new Vector[n];
	for(int i = 0 ; i < n ; i++)
	{
		a[i] = Vector(i+1);
	}
}

///@brief Constructor of a specific tensor
///@param n dimension of the created SymTensor2
///@param c type of SymTensor2 : k for kronecker
SymTensor2::SymTensor2(int n, char c)
{
	dim = n;
	a = new Vector[n];
	for(int i = 0 ; i < n ; i++)
	{
		a[i] = Vector(i+1);
	}
	switch(c)
	{
		case 'k':
		{
			for(int i = 0 ; i < n ; i++)
				for(int j = 0 ; j <= i ; j++)
				{
					if(i == j)
						a[i][j+1] = 1.;
					else
						a[i][j+1] = 0.;
				}
		}break;
	}
}

///@brief Copy constructor
///@param t SymTensor2 object to copy
SymTensor2::SymTensor2(const SymTensor2& t)
{
	dim = t.dim;
	a = new Vector[dim];
	for(unsigned int i = 0 ; i < dim ; i++)
	{
		a[i] = Vector(i+1);
		for(unsigned int j = 1 ; j <= i+1 ; j++)
		(a[i])[j] = (t.a[i])[j];
	}
}

///@brief Copy constructor for temporary objects
///@param t SymTensor2 object to copy
SymTensor2::SymTensor2(SymTensor2* t)
{
	dim = (*t).dim;
	a = new Vector[dim];
	for(unsigned int i = 0 ; i < dim ; i++)
	{
		a[i] = Vector(i+1);
		for(unsigned int j = 1 ; j <= i+1 ; j++)
		(a[i])[j] = ((*t).a[i])[j];
	}
}

///@brief Destructor
SymTensor2::~SymTensor2()
{
  if(dim != 0)
      delete[] a;
}

///@brief Overload of () operator
///@param i index of line
///@param j index of column
double& SymTensor2::operator()(unsigned int i, unsigned int j)
{
	assert((i <= dim)&&(i > 0));
	assert((j <= dim)&&(j > 0));
	if(j > i)
	{
		int temp;
		temp = j;
		j = i;
		i = temp;
	}
	return (a[i-1])[j];
}

///@brief Overload of = operator
///@param t SymTensor2 object to copy
SymTensor2& SymTensor2::operator = (const SymTensor2& t)
{
	assert(t.dim > 0);
	dim = t.dim;
	a = new Vector[dim];
	for(unsigned int i = 0 ; i < dim ; i++)
	{
		a[i] = Vector(i+1);
		for(unsigned int j = 1 ; j <= i+1 ; j++)
		(a[i])[j] = (t.a[i])[j];
	}
	return (*this);
}

///@brief Overload of = operator for temporary objects
///@param t pointer to SymTensor2 object to copy
SymTensor2& SymTensor2::operator = (SymTensor2* t)
{
	assert((*t).dim > 0);
	dim = (*t).dim;
	a = new Vector[dim];
	for(unsigned int i = 0 ; i < dim ; i++)
	{
		a[i] = Vector(i+1);
		for(unsigned int j = 1 ; j <= i+1 ; j++)
		(a[i])[j] = ((*t).a[i])[j];
	}
	return (*this);
}

///@brief Copies SymTensor2 without size check (use with caution)
///@param t SymTensor2 object to copy
SymTensor2& SymTensor2::operator &= (const SymTensor2& t)
{
	if(!t)
	{
		if(dim)
			delete[] a;
	}
	else
	{
		if(dim)
			delete[] a;
		a = new Vector[dim = t.dim];
		for(unsigned int i = 0 ; i < t.dim; i++)
			(*this).a[i] &= t.a[i];
	}
	return *this;
}

///@brief Cast operator to Tensor2
SymTensor2::operator Tensor2()
{
	assert(dim > 0);
	Tensor2 res(dim);
	for(unsigned int i = 1 ; i <= dim ; i++)
		for(unsigned int j = 1 ; j <= i ; j++)
		{
			if(i==j)
				res(i,j) = a[i-1][j];
			else
			{
				res(i,j) = a[i-1][j];
				res(j,i) = a[i-1][j];
			}
		}
	return res;
}

///@brief Return dimension
int SymTensor2::Dim()
{
	return dim;
}

///@brief Return deviator part of tensor
SymTensor2 SymTensor2::Dev()
{
	assert(dim != 0);
	double trace = 0.;
	for(unsigned int i = 1 ; i <= dim ; i++)
		trace += (a[i-1])[i];
	trace /= dim;
	SymTensor2 res(dim);
	for(unsigned int i =1 ; i <= dim ; i++)
		for(unsigned int j = 1 ; j <= i ; j++)
		{
			if(i ==j)
				(res.a[i-1])[j] = (a[i-1])[j] - trace;
			else
				(res.a[i-1])[j] = (a[i-1])[j];
		}
	return res;
}

///@brief Return trace of tensor
double SymTensor2::Trace()
{
	assert(dim != 0);
	double trace = 0.;
	for(unsigned int i = 1 ; i <= dim ; i++)
		trace += (a[i-1])[i];
	return trace;
}

///@brief Overload of * operator : product scalar-SymTensor2
///@param s left operand
///@param t rigt operand
SymTensor2 operator * (const double& s, const SymTensor2& t)
{
	assert(t.dim > 0);
	SymTensor2 res(t.dim);
	for(unsigned int i = 1 ; i <= t.dim ; i++)
		res.a[i-1] = s*t.a[i-1];
	return res;
}
///@brief Overload of * operator : product SymTensor2-scalar
///@param t left operand
///@param s rigt operand
SymTensor2 operator * (const SymTensor2& t, const double& s)
{
	return s*t;
}

///@brief Overload of + operator : addition of 2 Symtensor2
///@param t1 left operand
///@param t2 right operand
SymTensor2 operator + (const SymTensor2& t1, const SymTensor2& t2)
{
	assert(t1.dim > 0);
	assert(t1.dim == t2.dim);
	SymTensor2 res(t1.dim);
	for(unsigned int i = 1 ; i <= t1.dim ; i++)
		res.a[i-1] = t1.a[i-1] + t2.a[i-1];
	return res;
}
///@brief Overload of - operator : substraction of 2 SymTensor2
///@param t1 left operand
///@param t2 right operand
SymTensor2 operator - (const SymTensor2& t1, const SymTensor2& t2)
{
	assert(t1.dim > 0);
	assert(t1.dim == t2.dim);
	SymTensor2 res(t1.dim);
	for(unsigned int i = 1 ; i <= t1.dim ; i++)
		res.a[i-1] = t1.a[i-1] - t2.a[i-1];
	return res;
}

///@brief Overload of | operator : contracted product A_ij = B_ik C_kj
///@param t1 left operand
///@param t2 right operand
Tensor2 operator | (const SymTensor2& t1, const SymTensor2& t2)
{
	assert(t1.dim == t2.dim);
	Tensor2 res(t1.dim);

	for(unsigned int i = 1 ; i <= t1.dim ; i++)
		for(unsigned int j = 1 ; j <= t1.dim ; j++)
		{
			res(i,j) = 0.;
			for(unsigned int k = 1 ; k <= t1.dim ; k++)
			{
				if((i >= k)&&(k >= j))
					res(i,j) += (t1.a[i-1])[k] * (t2.a[k-1])[j];
				else if((i >= k)&&(k < j))
					res(i,j) += (t1.a[i-1])[k] * (t2.a[j-1])[k];
				else if((i < k)&&(k >= j))
					res(i,j) += (t1.a[k-1])[i] * (t2.a[k-1])[j];
				else
					res(i,j) += (t1.a[k-1])[i] * (t2.a[j-1])[k];
			}
		}

	return res;
}

///@brief Overload of || operator : double contracted product A = B_ij C_ji
///@param t1 left operand
///@param t2 right operand
double operator || (const SymTensor2& t1, const SymTensor2& t2)
{
	assert(t1.dim == t2.dim);
	double res = 0.;
	for(unsigned int i = 1 ; i <= t1.dim ; i++)
		for(unsigned int j = 1 ; j <= i ; j++)
		{
			if(i==j)
				res += (t1.a[i-1])[j] * (t2.a[i-1])[j];
			else
				res += 2. * (t1.a[i-1])[j] * (t2.a[i-1])[j];
		}
	return res;
}

///@brief Return true if SymTensor2 exists
///@param t SymTensor2 to be tested
bool operator ! (const SymTensor2& t)
{
	return !t.dim;
}

///@brief Overload of << operator
///@param s out stream
///@param t SymTensor2 object to put in the out stream
std::ostream& operator<<(std::ostream & s, const SymTensor2& t)
{
	for(unsigned int i = 0 ; i < t.dim ; i++)
	{
		s << t.a[i];
		if(i != t.dim - 1)
			s << "\n";
	}
	return s;
}


